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We study the dynamics of the mean field model of a Bose- 
Einstein condensed atom cloud in a parametrically forced trap 
by using analytical and numerical techniques. The dynamics 
is related to a classical Mathieu oscillator in a singular po- 
tential. It is found that there are wide resonances which can 
strongly affect the dynamics even when dissipation is present. 
Different geometries of the forcing are discussed as well as the 
implications of our results. 



A,(i) = Aj,o(l + eiCos(uJit)) 
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I. INTRODUCTION 

The recent experimental realization of Bose-Einstein 
condensation (BEC) in ultra-cold atomic gases, |l],|| has 
triggered the theoretical exploration of the properties of 
Bose gases. Specifically there has been a great interest 
in the development of applications which make use of 
the properties of this new state of matter. Perhaps, the 
recent development of the so-called atom laser || is the 
best example of the interest of these applications. 

The current model used to describe a system with 
a fixed mean number N of weakly interacting bosons, 
trapped in a parabolic potential V{f) is the following 
Nonlinear Schrodinger Equation (NLSE) (which in this 
context is called the Gross-Pitaevskii equation (GPE)) 



dtp 



h 

—V 2 ^ + V(r,t)^ + U \i;\ 2 i;, (1) 

which is valid when the particle density and temperature 
of the condensate are small enough. Here Uq = AttH a/m 
characterizes the interaction and is defined in terms of the 
ground state scattering length a. The normalization for 
ip is N — J \t/j\ 2 d 3 f, and the trapping potential is given 
by 

V(r, t) = \mv 2 (\l(t)x 2 + Xl(t)y 2 + A*(i)z 2 ) , (2) 



v, ( , (r) = x,y,z) are, as usual^functions that de- 
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scribe the anisotropics of the trap M. In real experi- 
ments with stationary systems they are constants and 
the geometry of the trap imposes usually the condition 
A x = A y = 1. A z = v z /v is the quotient between the fre- 
quency along the z-direction v z and the radial one v r = v. 

At this point we want to emphasize here that our anal- 
ysis does not restrict to the stationary case. Instead it fo- 
cuses on situations where A^ are periodic sinusoidal func- 
tions of the time. To be more precise we will adopt the 
notation 



with i = x, y, z. 

In its original derivation Eq. (]]]) was supposed to be 
strictly valid in the T = and low density limits. Recent 
theoretical work extends the applicability of the GPE to 
the high density limit . On the other hand linearized 
stability analysis based on perturbative expansions on 
\/N seems to point that the validity of the equation is 
restricted to the cases where no exponential separation of 
nearby orbits appear as it happens for example in chaotic 
pulsations of the atom cloud |7j|^] . 

A lot of work has concentrated in the analysis of the 
resonance structure when a periodic time dependent per- 
turbation is applied to the magnetic field, mainly be- 
cause of the availability of experimental data ||. The 
theoretical advances include both semi-classical analysis 
of Eq. (p]in the hydrodynamic limit |l0| , variational 
methods [11 12 and numerical simulations j^Jy|]. Re- 
cent research based on Eq. ([!]) also includes the study of 
the different problems using tools from nonlinear science 
|17 -[19| and the analysis of its multicomponent extensions 

[2g,|y. 

The variational approach provides us with very simple 
equations for the evolution of the widths of the Gaussian 
atom cloud together with a simple picture of the move- 
ment of the condensate p2| . However, those equations 
were initially derived for the static field case and only 
applied to the analysis of the normal modes and frequen- 
cies of the condensate motion in the weak perturbation 
regime. Similar equations are found using various scaling 
arguments |l^Jl^,0| or moment theory pCfl , 

Evolution of the condensate in a time dependent trap 
has been addressed in many different ways. There are 
papers where the GPE is solved numerically, either by 
watching the time evolution of a suitable initial condition 
1 13,[l4| or by using a linear expansion in normal modes 
\14\. And there are other works [p"5[|7[ where the au- 
thors derive a set of ordinary differential equations using 
scaling arguments plus the Thomas-Fermi approxima- 
tion. Among all the papers that treat time dependent 
potentials, the one which is closer to what we present 
here is 0. In that work the authors study the excita- 
tions of condensed and non-condensed clouds under a 
periodic parametric perturbation and find that the con- 
densate depletes and that chaotic evo lutio n comes out. 
We will comment more on this at Sec. 
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It is our intention in this paper to explore the behav- 
ior of the condensate under periodic perturbations of the 
trap strength using analytical and numerical tools. Ana- 
lytical techniques include exact results and a variational 
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analysis extended to situations with a time dependent 
potentials. Comparison of the results of both methods 
will allow us to state rigorously and derive a simple model 
for the resonant behavior and to qualitatively predict the 
evolution of the system under other conditions (less sym- 
metry dissipation regimes, non-condensed corrections, 
etc). Other conclusions as well as experimental implica- 
tions will be found in the course of the analysis. 
Our plan is as follows: In Sec. || 

we obtain a set of 

exact equations for the evolution of the center of mass of 
the condensate. The result is found to be three uncou- 
pled Mathieu equations. We show that a complete per- 
turbative analysis is possible and obtain the resonance 
structure for the center of mass. In Section [II we intro- 
duce the time dependent variatio nal model for the time 
dependent GP equation 



IV we solve the GPE 
We show the 
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numerically in the radially symmetric case 
resonance structure for the condensate in all regimes - 
small and large amplitude oscillations -. Next we turn 
to the variational equations for the widths and show that 
numerical simulations agree qualitatively with the pre- 
ceding picture. Finally, we use those ordinary differen- 
tial equations to predict the locations of the resonances 
in a simple model that connects the linear and the non- 
linear cases using tools from Sec. |l[ This is the main 
result of our paper. In Sec. we consider the possibility 
of extending the analysis to the multidimensional case 
and study the effect of the coupling between the differ- 
ent variational parameters. The resulting predictions are 
found to agree with a posterior set of 3D simulations of 
the GPE equation. In Sec. VI we comment on the effect 
of losses. We show that these resonances are of a per- 
sistent nature and show how chaotic evol ution may arise 
from losses effects. Finally in Sec. VII we discuss the 
experimental implications of our results and summarize 
our conclusions. 

Note that, unless otherwise stated, all magnitudes are 
in adimensional units. To adimensionalize these figures 
we have employed the change of variables that is intro- 
duced in Sec. 
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II. EXACT ANALYTICAL RESULTS 

A. Variational form of the GP equation 

It can be proved that every solution of Eq. ([!]) is a 
stationary point of an action corresponding, up to a di- 
vergence, to the Lagrangian density 

in f ar &i>\ 



Zm 



(4) 



where the asterisk denotes complex conjugation. That 
is, instead of working with the NLSE we can treat the 
action, 



S = J Cd 3 rdt = J 1 L(t)dt 



(5) 



and study its invariance properties and extrema, which 
are in turn solutions of Eq. ([[]). 

For instance, from the invariance of Eq. (|J) under 
global phase transformations, one can assure the conser- 
vation of the norm of the wave function 



N= / |V>r<rr, 



(6) 



which in this context is interpreted as the number of par- 
ticles in the Bose condensed state. We can also define 
another quantity, 



E = 



2m 



|W>| 2 + V(r, t)i> + 



d 3 r, (7) 



that can be thought of as the energy, and whose time 
evolution is simply 



— = / — \ihrd 6 r. 
dt / dt m 



(8) 



Thus, when the potential is not time dependent, the en- 
ergy is another conserved quantity. And when the po- 
tential has the form (Q) , the evolution of energy can be 
easily connected to that of the mean square radii of the 
cloud 



dE 1 

= —mv 

dt 2 ^ dt 

rf=x,y 7 z 



(9) 



Eqs. (g) and (g) are also useful to test the stability of 
the numerical scheme we will use to simulate Eq. ([!]). 



B. Newton's equations for the center of mass in a 
general GPE 

Let us consider the following function 

${r , )=<K?-?o), (10) 

where is a solution of Eq. |]. Substituting it in Eq. (||) 
and calculating the averaged Lagrangian, we obtain 



L = J Cd 3 r = Lf ree [4>] + L cm [(f>,r ]. 



(11) 



The Lagrangian has been splitted in two parts: one, 
the "free" contribution L/ ree , which depends only on <f> 



(12) 



Lfree[<f>] = J y (^t 0* " <f>* d t cj>) d\ 



and another one, L cmi that includes both the potential 
and the displacement r*o 
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,r = 



V{r + r n )\4>{r)\ 2 - ift^V^b} 



cfr. 



(13) 



If we impose that the action be stationary for some 
fo(t), i.e. if we use Lagrange's equations (p7|), we get 

|<- i W> = -<|-F(r + fo)>. (14) 

Here the brackets denote, as usual, the mean value of an 
operator over the unperturbed wave function, 0, < A >= 
/ <j)*(r)A(j)(r)d 3 r. 

As a final step let us use that is a solution of the 
GPE. Then Eq. (O) must be satisfied at least for r = 



d d 
— < -ihV > = -< ^—V{r) > 
at ori 



(15) 



It is now easy to prove the relation between the mean 
value of the moment operator, iftSJ , and the speed of the 
center of mass. We start from 



dt <r>= 



(16) 



Next we replace the time derivatives with spatial ones 
using Eq. (|l|) and its complex conjugate 
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dt ih I 2m 



V 2 4>-<t>—V 2 4>*Wr, (17) 



x + (1 + ecos(uit))x = 0. 



(20) 



This equation is known as Mathieu's equation. It is a 
well know problem which appears frequently in the study 
of parametrically forced oscillators and where one can 
obtain a lot of information by analytical means (2^ ||(]] . 

First, Floquet's theory for linear ODE with periodic 
time dependent coefficients shows that Eq. fl2C| ) 
has an infinite set of instability regions in the param- 
eter space. The limits of these zones can be found 
and have the shape of wedges that start on the points 
(w mi „,e min ) = (2,0), (1,0), (2/3,0), and widen as e is 
increased up from zero. Inside this regions at least one 
branch oscillates with an exponentially increasing ampli- 
tude. 

Either with an asymptotic method, or by making use of 
the singular perturbation theory, we can also locate those 
resonances and study the evolution of the system around 
them. For a perturbation frequency close enough to the 
first resonance, that is for \w — 2| = o(l), an asymptotic 
method j23j yields up to first order 



± 



4^2 



<5 2 , 



r~ce at cos(ujt/2 + 9 ). 



(21a) 
(21b) 



Here we see that for some values of 5 and e the exponent 
a is a positive real number and the amplitude of the 
oscillations grows unlimitedly. Also, the strength of the 
resonance is maximum for a value of 



and finally integrate this expression to obtain 



— < 7">=< > 
dt 



(18) 



Eqs. ( |15[ ) and (|18|) are the quantum equivalent of New- 
ton's Second Law and are exact for functions <f> that sat- 
isfy the GPE, in fact these equations coincide with the 
ones appearing in the linear Schrodinger equation. 



C. The condensate in a harmonic trap. Mathieu 
equations for the center of mass 

In our setup V(r) is a harmonic potential with trap 
strengths of the form of Eq. (|j|). Thus, Eqs. (|l5|) become 
a set of three decoupled ODE 



d 2 1 2,2, s 

< x >= --mw \ x (t) <x>, 
d 2 1 

<y>=-2 TOtj2A 'W <y >> 

d 2 1 

—t < z >= --mu 2 \ 2 (t) < z > . 
dt 2 2 



(19a) 
(19b) 
(19c) 



After a change of scale, all of the preceding equations are 
equivalent to a model one that we will write as 



-1 



e 2 



e 2 + 0(e 4 ). 



(22) 



A second order Taylor expansion in Eq. ( pia| ) lays the 
following limits 



\ u -2\ < - + —. 
1 1 - 2 32 



(23) 



The treatment of other resonances is more difficult as 
they are caused by higher order terms -at least of second 
order in the u> = 1 case-. In practice this means that 
they have a smaller region of influence and that they are 
not so strong. One has to choose large values of e, initial 
conditions (x, x) not too close to the equilibrium point 
and an excellent numerical integration method, in order 
to find real instabilities. However, if one concentrates 
not just in looking for exponential divergences, but on 
efficient pumping, it will be easily observed that on top 
of these subharmonics there are peaks of the energy gain 
speed (See Fig |b). 

Finally, we wish to point out that these resonances 
are of a peculiarly persistent nature: as we will precise 
in section VI, they do resist even the presence of dis- 
sipation. This has a serious and inmediate consequence 
which is that feeding the condensate in a resonant regime 
can result in a large amplitude oscillation of the center of 
mass. On the other hand, as we will outline later in Sec. 



3 



VI, a measure of this effect can give us some insight in 
the losses effects, as well as of any additional terms that 
could be added to Eq. (|j) so as to better modeling the 
condensate. 



III. VARIATIONAL EQUATIONS FOR THE 
TIME-DEPENDENT GP EQUATION WITH A 
PERIODIC PARAMETRIC PERTURBATION 

Although the center of mass of the wave packet satis- 
fies very simple equations it does not happen the same 
to other parameters such as the width. Only in the two 
dimensional case it is possible to apply moment tech- 
niques and find analytically its time evolution |^7],^8| but 
in the fully three dimensional problem no exact results 
have been derived yet using the moment technique. 

To simplify the problem, we restrict the shape of the 
function ip to a convenient family of trial functions and 
study the time evolution of the parameters that define 
that family. A natural choice, which corresponds to the 
exact solution in the linear limit (Uo = 0) and provided 
quite good results in our previous works lllUT2] is a three 
dimensional Gaussian-like function with sixteen free pa- 
rameters 



tp(x,y,z,t) 



n 



A | | exp 

77= x,y,z 



-[v - vo? 

2w 2 



(24) 



For a matter of convenience and ease of interpretation 
we will make here a change of parameters, from A and A* 
to N (the norm of the wave function) and <j> (its global 
phase): 



A 



N 



n 3 / 2 w x w y w z 



J4> 



(25) 



The rest of the parameters are w v (width), a„ (slope), 
(square root of the curvature radius) and 770 (center 
of the cloud) . 

This trial function must now be placed in Eq. (^|) to 
obtain an averaged Lagrangian per particle 



Cd 6 r 



+ E 1 -f + % 2 \ \ h ^ + it ft + ^ 2A ' w 



L _ 1 
N ~ N 

2h 2 

m 



-y\— 

m ' 1 2w 2 



Uq 



N 



W x WyW z 



(26) 



The evolution of the parameters is ruled by the corre- 
sponding set of Lagrange equations 



d ( dL 
dt \dqj 



dL_ 



(27) 



which give us equations for the conservation of the norm, 
dN 



dt 



= 0; 



the movement of the center of mass, 
770 + mv 2 X n (t)r] = 0; 
the evolution of slope and curvature, 

a mtbr, 



2hw v ' 



a,, = — Vo - ZPvVo; 
and finally the evolution of the widths, 



Wy + V 2 \ 2 (t)Wy = 



,2 \2/ 



h 2 


1 


U N 


m 2 


w l 


2\/2m TT 3 / 2 w 2 w y w z 


h 2 


1 


Uo N 


m 2 


w l 


2y / 2m tt 3 / 2 w x w 2 w z 


h 2 


1 


Uo N 


m 2 


wl 


2\/2m n 3 / 2 w x w y w 2 



(28) 
(29) 

(30a) 
(30b) 

(31a) 
(31b) 
(31c) 



As one can see, the introduction of a time dependent 
potential does not affect the form of the equations, which 
remain the same as those of jl2| . 

Let us introduce the constants P = \f2/irNa/ao 
(strength of the atom-atom interaction) and ao = 
y/h/ (mv) (harmonic potential length scale), as well as a 
set of rescaled variables for time, r = 1st, and the widths, 



io v rj, (v = x i V, z )- This leads us to 



v x + X 2 x (t)v x = \ + 



P 



V 3 V 2 VyV Z ' 



X 2 y (t)Vy 



P 



V 3 V X VyV 2 



(32a) 
(32b) 
(32c) 



This set of rescaled variables and units is the one that we 
will use throughout the paper, unless otherwise stated. 



IV. ANALYSIS OF THE RADIALLY 
SYMMETRIC CASE 

A. The equations 

In this section we will analyse the case where the trap 
has the same strength on all directions, that is 

V(x, y, z) = l -mu 2 \ 2 (t)(x 2 + y 2 + z 2 ). (33) 
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and the solutions are supposed to be radially symmetric. 
This high degree of symmetry simplifies the equations 
considerably. First, the variational model for the widths 
of the condensate (31a) reduces to a single ODE for the 
radial width v(t) 



1 



P 



v = -\ 2 (t)v + ^ + -. (34) 

And secondly the following change of function 

u(r) 



ip(r,t) = A- 



(35) 



with the constrains 

u(r) — ► 0, r — > 
u{r)dr = 1 



(36a) 
(36b) 

\A\ 2 = AttN, (36c) 
transforms Eq. (uh into the one-dimensional PDE 



ihd t u = -—dlu + <^ -mv 2 \ 2 (t)r 2 + 4tt-^— L 4- 
2m 2 47T H 



(37) 



This is the equation that we have actually solved numer- 
ically. 



B. Numerical study of the equations 



The numerical solution of Eqs. ( |34| ) and (37) is not 
a trivial job. We can see without much effort that both 
equations are stiff Q when the width of the cloud be- 
comes very small, due to the presence of strong singular 
potentials. Thus, we need numerical methods that ac- 
count for the nonlinearities and are stable enough to be 
trusted when close to a resonance. 

To solve Eq. ( j34j) we have used an adaptive step 
size Runge-Kutta-Fehlberg method, a Dormand-Princc 
pair ||, the ODE Suite of MATLAB @, and finally 
Vazquez's conservative scheme fl3l| ] -a finite differences 
scheme that conserves a discretized version of the energy 
and is unconditionally stable-. All of them gave the same 
accurate results about the frequencies and amplitude of 
the oscillations, the regions of divergence, etc. 

To solve Eq. (37) we have utilized a modification of a 
second order accurate finite difference scheme developed 
in ]33| . This new scheme is time reversible, conserves the 
norm and has a discrete analogue for Eq. (||) which pro- 
vides enhanced stability |Q. On the other hand, even 
using the best methods, we face another important diffi- 
culty which is the finite size of either the spatial grid -in 
finite differences schemes- or the momentum space -in 
spectral or pseudospectral methods-. This size effect be- 
comes particularly important in the case of parametrical 



perturbations and imposes a severe limit on the time for 
which simulations may be trusted. 

Using all this computational machinery we have 
achieved several important results. First we have checked 
our programs with low amplitude oscillations. In the 
PDE we imposed a gaussian initial condition of width vO 
and used this same value as initial condition for Eq. (|34|) . 
By this procedure we obtained the linearized excitation 
frequencies of the condensate, concluding that there is 
a significant coincidence of both models as was shown 
in p^ |. In the large amplitude simulations we found 
a surprising result which is that the variational model 
still follows the evolution of the PDE with a 90% of ac- 
curacy even in situations where the condensate remains 
no longer gaussian but gains an important contribution 
from the first and second modes. Another confirmed fact 
is that the width of the condensate develops fast tough 
bounces against the "origin" (Fig. 0a-b), with its excita- 
tion frequency deviating from the linearized predictions 
of |l2[ | and approaching those of the harmonic trap. 

We have also simulated the system with a periodic time 
dependent perturbation of the form 



\ 2 (t) = 1 + ecos(ojt). 



(38) 



Using the variational equation Eq. |34| we have scanned 
the parameter space (w,e). We have found one region 
where the radial width diverges exponentially (Fig. |2ja) 
and causes most numerical methods to fail after a finite 
time, and two more where the width grows almost ex- 
ponentially up to a point where the simulation cannot 
account for this growth. All of these zones have the form 
of wedges, with a peak at (w m i n ,e m j„), and a growing 
width as e is increased. The most important resonance 
stands on u! m in = 2.04 (See Fig. |||b). The other ones are 
weaker and rest on w m , n = 1.02,0.68. We have studied 
a wide range of setups and found that these frequencies 
change no more than a 0.5% depending on the initial 
conditions and the nonlinearity. On the other hand, the 
lowest perturbation amplitude for which the resonance 
exists t m in does exhibit a strong dependence on the ini- 
tial conditions. Indeed, around the weaker resonances, 
instability is never reached for points close to the mini- 
mum of the potential. However, as we already pointed in 
Sec. H this is probably a numerical effect. 

In order to study the location of the resonances we have 
also made some plots of the efficiency of the energy ab- 
sorption process against the perturbation frequency for 
different values of the perturbation amplitude for the 
variational system and the partial differential equation 
(Figs. |]a-b, ||a-b). The way we have measured "effi- 
ciency" is by letting the system evolve for a fixed time 
and then computing the maximum value of the mean 
square radius of the cloud. For the sake of simplicity and 
to approach the experimental setups, we have used an 
equilibrium state of the static GPE as initial condition. 

In a rather complete inspection of the parameter space 
using the efficiency plots, we have found that, though the 
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resonance regions cannot be precisely delimited because 
of the dependence on the initial data, they do exist and 
behave much like the variational model predicted. As we 
see m Figs. | and | there are two important features 
in the response of the condensate. The main one is the 
width of the resonances and the dependence of that width 
on the strength of the perturbation. The second impor- 
tant feature is the change of the frequency for which the 
perturbation is most efficient. This peak is centered on 
the frequency of the linearized model only for very small 
perturbation amplitudes, and switches to the trap nat- 
ural frequencies very quickly as the amplitude is made 
stronger - of about 10% or so-. For even stronger ampli- 
tudes the optimal frequency decreases slowly. 

Another consequence of this work is that response of 
the cloud is stronger in the PDE than in the variational 
simplification. For instance, it is possible to find (See 
Fig. ||a) perturbation amplitudes that do not cause an 
significant growth in Eq. (34|) but make the cloud width 
increase quite linearly in the PDE. We have also stud- 
ied the case where an exponential growth of the width is 
present in both equations -([Mj) and (|37])- and we have 
seen that the growth is qualitatively similar though there 
is a tendency in the exact model to exhibit slightly larger 
amplitude motion. These discrepancies are originated by 
high modes which are not present in the variational treat- 
ment. In fact, even though an important part of the cloud 
remains close to the origin, it is observed that for long 
times a long tail appears which is hard to appreciate but 
is responsible for the growth of the mean square width 
(and the presence of higher order modes) . 



C. Analysis of the Mathieu equation with a singular 
potential 

In this subsection we are to develop a simplified model 
that explains why do resonances appear in the perturbed 
GPE equation. This model makes use of the variational 
equations for the cloud width but does not care for the 
actual shape of the variational ansatz. The reason for it 
is that both the approximate model and the exact one 
agree for short times, the response of the latter being 
always stronger. 

Let us limit ourselves to Eq. (|34|). In the previous 
subsection we said that this system is very stiff and that 
the origin acts as an elastic wall. In view of this, it is 
intuitively appealing to replace the singular but differen- 
tiable potential 1 /v 3 + P/ v 4 with a discontinuous bounce 
condition on the origin, i.e. an impact oscillator. Nu- 
merical simulations confirm that this approximation is 
good for large amplitude oscillations. Also, the connec- 
tion between soft singular potentials and impact oscil- 
lators has been widely studied and these kind of models 
have found ample application in real-life situations in the 
fields of Physics and Engineering (See p5| , p6| for many 
references) . 



Since the variational model is a lossless one, our bound- 
ary condition must be elastic. We replace Eq. (34) with 
the following one 



lim (v, v) 



+ X 2 (t)v 
(0 + ,V c ) 



o, 

lim (v, v) 

t~*tt 



(39) 



(0 H 



-v c . 



where t c denotes any isolated instance when the system 
bounces against the v — singularity. 

Let us show that this equation is in turn equivalent 
to an elastic oscillator without barrier conditions. We 
introduce the change of variables 

v=\u\, (40) 

where u is an unrestricted real number and satisfies the 
following one-dimensional harmonic oscillator equation 



u+ \ 2 {t)u = 0, 



(41) 



It is now easy to prove that every solution of Eq. ( |4l| ) 
provides a solution of Eq. (|39|). And vice versa, from 
every solution of Eq. ( |39| ) it is possible to construct a 
solution of Eq. ( pd] ) , unique up to a sign. 

So, what do we have now? We have proved that for 
large amplitude motion Eq. (|34| ) behaves as Mathieu's 
equation. This implies that in the regime of medium 
to large amplitude oscillations, or in a situation of large 
amplitude perturbations, Eq. ( |34"| ) will have instability 
regions that are more or less centered on Mathieu's fre- 
quencies u> = 2, 1,1/3,.... This prediction is indeed con- 
firmed by the numerical simulations: the main resonance 
is capable of causing an exponential divergence, while 
the other ones are harder to track down but do appear 
in plots of pumping "efficiency" . 

Another, important but minor result of this equiva- 
lence is that these resonance regions must become wider 
and move on to smaller frequencies as the perturbation 
amplitude is increased. This result is also obtained in 
the numerical simulations. Indeed, from the graphs we 
can estimate how fast the optimal frequency decreases 
as e grows, and we will see that the order of magnitude 
corresponds to that of Eq. ( |22| ) . 

Summing up, what we have seen here is that for low 
amplitude oscillations the condensate moves in an effec- 
tive potential that is parabolic. Thus, the frequency that 
excites the condensate most efficiently in a parametric 
way is the one that results from the linearization of Eq. 
(p4|). On the other hand, when we start to consider large 
amplitude oscillations, we find that the harmonic trap 
gains importance over the details of the well. It is in this 
situation that the instability arises, and it happens for 
perturbations that oscillate according to multiples of the 
frequency of the trap. 

Even more, higher modes are little influenced by the 
nonlinearity, just because they are more spread and the 
value of \ip(r, t)\ 2 is smaller. As a consequence, the en- 
ergies of these modes come even closer to those of the 
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linear harmonic oscillator, resulting in the fact that the 
response of the GPE is stronger than what the variational 
model, limited to a gaussian shape, predicts. 



V. ANALYSIS OF THE NONSYMMETRIC CASE 

After finishing the study of the radially symmetric 
problem, it seems a natural step to proceed with the 
nonsymmetric one. However this step is not simple for 
several reasons, the main one being that the simulation 
of full three dimensional GPE is computationally very 
expensive work. Thus, it would be unwise to directly 
attack the full problem without gaining some insight on 
what is to be expected from the simulations by cheaper 
means. 

In our case the cheapest tool is the variational model. 
We have seen that it describes rather well the behaviour 
of a radially symmetric condensate. And, as we pointed 
out before, there are exact analytical studies [ p7|j28| 
where the moment method reduce exactly the two di- 
mensional GPE to a set of ODE which are similar to the 
ones we have. 



A. Predictions of the variational model 

When we remove the radial symmetry in the varia- 
tional equations, we are left with two to three coor- 
dinates, and the perturbation can bear many different 
forms. However, trying to follow the experimental setups 
p2| , we should once more take a sinusoidal time depen- 
dence for every X v (t) coefficient, as in Eq. (g). This 



A: 



())/ 



2,1,2/3,. 



(43) 



e yi e z = 0) from 



choice accounts both for the m = mode (e x 
0) and the m = 2 perturbations (e x 
the JILA experiment |2^] . In the latter case the potential 
is a parabolic one, with fixed frequencies on a rotating 
frame. However, for our trial function (Eq. (|lj) it behaves 
just as a potential of the form of Eq. (||)). 

Su bstit uting our effective perturbation frequencies into 
Eq. ( |32a| ) we get a set of three coupled Mathieu equations 
with a potential that is singular on the v x = 0, v y = and 
v z = planes. The singularities are at least as strong as 
l/v 3 , and the numerical simulations again confirm that 
they act as elastic walls, so we now proceed with a change 
of variables formally equivalent to that of Eq. (40): 



(42a) 
(42b) 



for r) = x, y, z. 

Now the situation is a bit more complex. The first 
new feature is the existence of several sets of instability 
regions. Due to having three a priori different constants 
\o v , the three oscillators in Eq. fll^ ) are not equivalent 
and we may get three sets of resonances in the (e^,w) 
space, each one containing the instability regions that 
start on the frequencies 



Numerical simulations of the variational equations for the 
m = type excitation confirm this prediction with a 
relative accuracy around 0.5% in the frequencies. The 
results show again that, opposite to the pure Mathieu 
equation, these "wedges" rest on a nonzero value of the 
perturbation amplitude e m m,7j- In Fig. we see the 
evolution of the condensate width with parameters close 
to the main resonance region. 

Another new feature is the possibility of coupling be- 
tween the widths of the condensate. This coupling is 
seen both in the m — and m — 2 excitations setups. In 
the first one, the perturbed width feeds the unperturbed 
one. Figure || demonstrates that efficiency is not very 
high. In the m = 2 case two widths are associated to 
the same trap frequencies and the perturbations are of 
equal magnitude and opposite sense. The explanation is 
that both widths block each other, eliminating the reso- 
nance and leaving just a bounded "movement" of small 
amplitude. 

As a side note, we must say that the variational model 
itself predicts the lack of resonances in the m = 2 setup. 
We already pointed out that the JILA m = 2 perturba- 
tion corresponds to a situation where the potential of the 
trap is rotated while maintaining its shape. It can be eas- 
ily proved that if we choose any family of trial functions 
with enough degrees of freedom for a general rotation, 
the variational solution will always stick to the poten- 
tial and rotate at the same speed. This also implies that 
the trial function (|24|) is not suitable for describing the 
condensate when this kind of perturbation is applied. 



B. Numerical study of the three dimensional GPE 

To perform the simulation of the full GPE we have used 
a fourier pseudospectral method, using typically a grid 
of 108 3 collocation points and integrating in time with a 
second order, symmetrized split step method p9| , [iT)[ . 

Our numerical study is based on a 0({/S.tf) scheme 
which behaves extremely well for long time runs. How- 
ever, as in the radial PDE, no matter how accurate the 
scheme is, there's a limit in the time during which sim- 
ulations can be trusted and this limit is imposed by the 
growth rate of the condensate width and the size of the 
grid. This limit is specially important for the m = per- 
turbation, where the condensate develops a large tail in 
the unperturbed direction, and this tail breaks the sim- 
ulation after a certain time. The more asymmetric the 
trap is, the sooner this effect comes out. In our simu- 
lations the grid was a box of 108 3 equally spaced points 
and whose sides measured from 20 to 40 length units, de- 
pending on the symmetry and intensity of the trap. This 
allowed us to track the condensate for about 12 periods in 
truly resonant setups, and for much more in nonresonant 
ones. 
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We have applied the algorithm to many different prob- 
lems. First we checked our programs against stationary 
and radially symmetric problems. Secondly we intro- 
duced time depe nden t traps and reproduced the calcu- 
lations of Section 



IV 



In both cases we got the expected 

results. 

The third set of experiments consisted in a resonant 
time dependent radially symmetric trap applied to sev- 
eral slightly asymmetric gaussian wave packets. The ini- 
tial asymmetry was slight enought to treat it as a weak 
perturbation and we saw that it departured little from 
the symmetric case -i.e., no modes with higher energy or 
angular moments break the exponential growth. 

The fourth and probably most important group of sim- 
ulations consisted in the study of the m = perturba- 
tion from the JILA p2| experiment. Here we confirmed 
the predictions of the previous subsection, that is we ob- 
tained at least one resonance region where the variational 
model predicts. 

We also observed that the response of the condensate 
as modeled by the GPE is stronger and exhibits a slightly 
more intense growth rate than what the variational model 
predicts. This and other similar results from the radial 



case (Sect. IV) favor the theory of a cooperative staircase 



effect, where the higher modes contribute to the energy 
absorption process without interfering in the evolution 
of the width. Figs. || and ^| show just two examples of 
the kind of evolution of have seen for this perturbation 
model on the resonant regions. 

The reason for this cooperative effect seems to lay in 
the energy level structure of the GPE equation. We 
have studied the evolution of the correlation a conden- 
sate wave against its initial data. In all cases the initial 
data was a displaced and deformed gaussian cloud, while 
the environment corresponded to a stationary trap. In 
the linear case it is easy to show that the spectrum of 
the correlation must reveal a subset of the eigenvalues of 
the hamiltonian, and indeed that is what we got. In a 
nonlinear context it is not clear what the frequencies of 
the correlation mean -they may be or may be not eigen- 
values of the GPE-, but at least we know that they must 
rule the energy absorption process somehow. What nu- 
merical experiments show is that for an extensive family 
of initial conditions these generalized "spectra" can be 
approximated by the formula E n — ujn + Eq, where E 
depends on the nonlinearity and the level spacing is a 
regular, harmonic one. 

This picture of equally spaced levels is intuitively ap- 
pealing to explain the existence of such strong reso- 
nances. On one side we have the fact that, in any other 
system, a continous parametric perturbation would be- 
come bounded as the higher state become populated. 
This is probably what one would first think when fac- 
ing this system. On the other side we see that due to 
being equally spaced these higher energy states are them- 
selves sensitive to the same resonant frequency and offer 
no resistence to the particle promotion process. In the 
end, this is more or less what happens in the variational 



equations, considered from the point of view of Classical 
Mechanics. 

However, it is not relevant for the existence of reso- 
nances whether they cause a sustained growth or not. 
What is confirmed without any doubt is that the pertur- 
bations that are most efficient in the variational model 
are also the most efficient in the full GPE. 



VI. ANALYSIS OF THE EFFECT OF LOSSES 

We now want to show the effect of a dissipative term in 
the variational equations. This term will be introduced 
in a phenomenological way so as to model the damp- 
ing of the oscillations of the condensate in regimes where 
the number of noncondensed particles is small. We will 
choose a viscous damping term that models well the be- 
haviour of the condensate in the experiments p2| . This 
choice introduces a significant loss of energy in the oscilla- 
tions while preserving the number of condensed particles. 
Using this term we will see that the resonance regions of 
the pure Mathieu persist. 

Let us see what Eq. (^0|) looks like once we add damp- 
ing: 



x + (1 + e cos(ujt))x + 72; = 0. 



(44) 



A simple change of variables x(t) — p(t)e~ l1 makes this 
new term disappear, transforming it back to a pure Math- 
ieu equation 



p + (1 - 7 2 + e cos(wi))p = 0. 



(45) 



With the introduction of the damping we are shifting the 
resonances to lower values, given by 



OJ 



2,1,2/3... 



(46) 



where v (7) = 1 — j 2 is the new effective frequency for 



the trap. We can approximately solve Eq. 
the first resonance, obtaining 



x(t) 



>-7)' COS (- 



"J) around 



(47) 



where A is given by Eq. (21a). 

This shows that the resonance regions in the parameter 
space are constrained to values of (uj, e) for which the 
strength of the resonance, A, is larger that the strength 
of the dissipative term. The new regions have a larger, 
nonzero, value of t m im and are typically thinner, but do 
not disappear unless 7 is very large. 

This effect is reproduced in the variational model when 
we introduce similar viscous damping terms. For in- 
stance, taking the data from the JILA experiment p2[ , 
we can estimate a condensate lifetime of about 110 ms 
and a value for 7 of 0.15 in natural units of the conden- 
sate. Such damping makes the e m i n value raise from 0.09 
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to 0.18 for the P = 9.2 case. Thus, the instability should 
not be appreciated unless the perturbation amplitude ex- 
ceeds the 20%. 

An interesting effect of damping is that the evolution of 
a continuously perturbed condensate outside the instabil- 
ity regions becomes more ordered than in the undamped 
mode since the motion is constrained to limit cycle syn- 
chronized to the frequency of the •parametric ■perturbation, 
and with a size that depends only on the perturbation pa- 
rameters, (u>, e). In Fig. |] we show the different limit cy- 
cles that appear under peridic perturbations. The largest 
one is always the one with its frequency on top of the 
peak of the resonance as shown by Fig. || and the size 
of the limit cycle decreases as the frequency is detuned 
from this value. 

Finally we wish to point out that the appearance of a 
limit cycle opens the door to a wide family of phenom- 
ena, from chaotic motion to bifurcation theory p6| , |37|] . 
This limit cycle would exist under a great variety of dis- 
sipative terms, and is not exclusive of linear damping. 
Also, the dependence of the limit cycle on the damping 
constant can be useful from the experimental viewpoint 
to separate the condensed and noncondensed clouds as 
will be pointed out in the last section. 



VII. CONCLUSION AND DISCUSSION. 

In this work we have analysed the resonant dynamics 
of the parametrically forced time dependent GPE using 
exact analytical techniques, approximate time dependent 
variational techniques and numerical simulations. All the 
results point to the existence of various resonant behav- 
iors associated to the same parameter regions concerning 
the motion of the center of mass and the width oscilla- 
tions of the wave packet. The role played by the nonlin- 
earity is to provide a strong repulsive term at the origin 
which acts as a barrier, which depending on the dimen- 
sionality and the symmetry of the external forcing could 
be stronger than the repulsive term related to the linear 
dispersion (kinetic energy term). 

We have developed a simplified version of the varia- 
tional equations for the spherically symmetric condensate 
under periodic sinusoidal change of the trapping poten- 
tial. This model is based on an impact oscillator, i.e. 
a harmonic oscillator with an elastic barrier condition 
which allows to find explicit solutions. When applied 
to the parametrically perturbed condensate, our model 
predicts that medium to large oscillations approach the 
harmonic trap frequencies, not the ones resulting from 
the linearization of the variational equations. The model 
also predicts that the response of the condensate to an 
external perturbation is ruled by the harmonic trap fre- 
quencies. In the case of a sinusoidal parametric pertur- 
bation it accurately predicts the existence of a family of 
resonances on multiples of the trap natural frequencies. 

In the end, what we get from this work is a precise pic- 



ture of the resonances for a wide family of equations that 
include the Eq. ([[]) and the harmonic oscillator. In this 
description for the radially symmetric case we have one 
base frequency that is essentially the same for both equa- 
tions (linear and nonlinear one) and which corresponds to 
the energy separation between the ground state and the 
first excited state of the harmonic oscillator up to a high 
precision. The invariance of this base frequency has been 
checked for a nonlinearity constant going from P = 9.2 
-the JILA |^2j experiment- to 20 times this value. It 
differs from the predictions of p8[ , where formulas re- 
garding the P — > and P — * oo are derived. Also there 
is a whole set of subharmonics of this frequency, all of 
which are capable of exciting the cloud quite efficiently. 
At least three of them have been found, both with signif- 
icant responses. These subharmonics are not predicted 
in ||. 

We have also found that for this kind of parametric 
drive the resonances are wide. The width grows with 
the strength of the interaction and decreases with the 
effect of dissipation. Both facts can be checked in the 
experiments by forcing the system for a longer time than 
what it is currently done. 

Both predictions are exact in the linear limit, Uq = 0, 
and have been confirmed with simulations coming both 
from the exact variational model and the GPE. The dis- 
cretization of the GPE has also allowed us to scan the 
parameters space, studying the efficiency of the pertur- 
bation process. In this study we have only found peaks 
centered on Mathieu's frequencies. 

For the general case with more than one degree of free- 
dom (axial and nonsymmetric cases) we obtain a set of 
two to three decoupled pure Mathieu equations. We have 
shown that, due to having more than one frequency, the 
predicted Mathieu resonances do exist in a larger num- 
ber. On the other hand, we have also seen that some 
of these resonances may disappear due to the locking of 
"equivalent" variables, an effect that our decoupled equa- 
tions do not account for. 

This resonance scheme for the nonsymmetric case has 
been confirmed with accurate simulations of the full GPE 
for the m = perturbation. An analysis of the correla- 
tion of a state against its initial data shows that both 
the linear and the nonlinear problems exhibit a spectral 
structure which is likely to present such behavior. 

Finally, damping has been shown to limit the effect of 
the parametric perturbation. Once more we have proved 
that only frequencies close to the Mathieu resonance re- 
gions do excite the condensate as a whole in an efficient 
way, causing the appearance of a stable limit cycle. All 
other frequencies are inefficient in the sense that the sys- 
tem stays extremely close to the equilibrium configura- 
tion, which acts as a focus. 

A main conclusion of this work is that for this set of 
resonances to exist one only needs a singularity that pre- 
vents from collapse. The variational method showed that 
the kinetic terms in the evolution equations guarantee a 
1/x 3 singularity as far as we impose a repulsive interac- 
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tion between the atoms in the cloud. This is the reason 
why we say that we have a family of systems that be- 
have much the same. An immediate result of this is that 
the response of the noncondcnsed atoms under the para- 
metrical perturbation will be qualitatively similar to that 
of the condensed ones, with the only difference that the 
former are subject to a more intense dissipation. But 
as we already saw in Sec. VI this dissipation can be 
enough to distinguish both kind of fluids: while the con- 
densed part might suffer an exponential growth, the un- 
condensed part might develop low amplitude bounded 
oscillations. 

We have also demonstrated that these resonances show 
up in the movement of the center of mass as well, causing 
any initial displacement of the center of mass to be ex- 
ponentially amplified while the perturbation works. Op- 
posite to our models for the widths, this is an exact pre- 
diction based solely on the GPE, and it shows that the 
parametrical perturbation may also have a disturbing ef- 
fect in the experiments. On the other hand, a measure 
of this effect can give us information about the intensity 
of dissipation and collision effects. 

All of the preceding statements are based solely on the 
GPE. In a few words, they include the existence of reso- 
nance regions both for the widths and the center of mass, 
the shape and the location of those regions, and its inten- 
sity as possible measure of damping. The failure of any 
prediction should be interpreted as a failure of the GPE 
to describe the condensate. Thus we provide with sim- 
ple experimental checks to perform a quantitative study 
of the regimes for which the GPE properly describes the 
Bose-Einstein condensates in time dependent traps. 

Finally, we must mention that throughout this work 
we have concentrated on regular motion regions in the 
parameter space. These regions can be "safely" reached 
in the experiments. There are many other cases where 
chaos appears in the variational equations and complex 
behavior is seen in the numerical simulations of Eq. (^) . 
Although the study of those disordered regions could be 
interesting from the nonlinear science point of view they 
seem not to be of interest for Bose-Einstein condensation 
since the exponential separation of nearby orbits which is 
characteristic of chaotic behavior has been shown to in- 
duce instabilities and take the system out of the regime 
where it can be described using the mean field GP equa- 
tion 
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FIG. 1. Radially symmetric condensate with P — 9.2 as 
described by the variational equations. Phase space picture 
for different small to large amplitude oscillations. 

FIG. 2. Radially symmetric condensate with P = 9.2 
and initial condition v = 1.6, subject to a periodic perturba- 
tion. Evolution of the width in the variational model (dashed 
lines) and in the GPE (solid line) when the perturbation is 
(a) (e = 0.15, oj = 4.00) and (b) (e = 0.2, oj = 4.00). (c) 
Instability region for the variational model in the parameter 
space around the main resonance. 

FIG. 3. Radially symmetric condensate with P = 9.2. 
Plot of the maximum amplitude of oscillations for the GPE 
after 40 time units. The initial condition corresponds to gaus- 
sian of width v = 1.6. Each line corresponds to a different 
value of e, from 0.05 to 0.3 in steps of 0.05. The frequency 
range covers (a) the main resonance and (b) the second im- 
portant one. 



FIG. 4. Radially symmetric condensate with P = 9.2. 
Plot of the maximum amplitude of oscillations for the vari- 
ational model after 40 time units. Initial conditions are 
v — 1.6, v — 0. Each line corresponds to a different value 
of e, from 0.05 to 0.3 in steps of 0.05. The frequency range 
covers (a) the main resonance and (b) the second important 
one. 

FIG. 5. Evolution of a cylindrically symmetric condensate 
(variational model, P = 9.2) under a sinusoidal perturbation 
{ui,e) = (2.04,0.1) of only the radial strength of the trap. 
Both (a) the radial v r and (b) the axial v z widths are plotted. 

FIG. 6. Phase space picture for the width of radially 
symmetric condensate subject to damping plus a periodic 
perturbation. The nonlinearity is P — 9.2, the damping 
7 = 0.15. The perturbation has in all cases e = 0.1, while 
the frequencies are, from the outer cycle to the inner one, 
lj = 2.15,2.4,3.0,4.0. 



FIG. 7. Evolution of a condensate in a cylindrically sym- 
metric trap (P — 9.2, A r = 1, X z = 2) subject to a m = 
perturbation (ui r ,e r ) — (2.00,0.15). Both the radial width v r 
(solid line) and the axial v z width (dashed line) are plotted. 

FIG. 8. Evolution of a condensate in a spherically sym- 
metric trap (P = 9.2) subject to a m = cylindrical pertur- 
bation (uj r , e r ) = (2.00, 0.15). Both the radial width v r (solid 
linea) and the axial width v z (dashed line) are plotted. 
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